Critical Step Length as an Indicator of Surface Supersaturation during Crystal Growth from Solution

The surface processes that control crystal growth from solution can be probed in real-time by in situ microscopy. However, when mass transport (partly) limits growth, the interfacial solution conditions are difficult to determine, precluding quantitative measurement. Here, we demonstrate the use of a thermodynamic feature of crystal surfaces—the critical step length—to convey the local supersaturation, allowing the surface-controlled kinetics to be obtained. Applying this method to atomic force microscopy measurements of calcite, which are shown to fall within the regime of mixed surface/transport control, unites calcite step velocities with the Kossel–Stranski model, resolves disparities between growth rates measured under different mass transport conditions, and reveals why the Gibbs–Thomson effect in calcite departs from classical theory. Our approach expands the scope of in situ microscopy by decoupling quantitative measurement from the influence of mass transport.


Finite Element Method
The relationship between the flow rate and the concentration boundary layer established at the surface of a crystal during a typical in situ AFM experiment was modelled using the finite element method in COMSOL Multiphysics (v5.6), applying the Transport of Diluted Species and Laminar Flow modules.

Geometry
The AFM experiments S1 employed a Bruker Nanoscope AFM that features a Digital Instruments flow cell. For the finite element model, we used the Digital Instruments MTFML flow cell geometry from a previous study, S2 with the addition of a 2 × 2 × 0.5 mm 3 rhombohedron S2 at the origin to represent the calcite crystal. The cantilever was based on the geometry of the probe described in ref. S1 and was positioned at the surface of the rhombohedron. A 9 × 9 µm 2 square patch on the crystal surface, with the AFM tip at the centre, was taken to represent the scan area. To accurately simulate the conditions near the cantilever, a multi-scale approach was adopted. A fine mesh of size 500 nm was used within a small domain encompassing the rhombohedral crystal and the cantilever, and a coarser mesh applied elsewhere. The meshes were refined until the solute flux within the region of interest changed by less than 3%.

Speciation Model
The model describes the transport of Ca 2+ and CO 2− 3 ions only, where the inlet concentrations were based on the highest-supersaturation experiment in ref. S1 In preliminary work, we examined the role of other species (HCO − 3 , NaCO − 3 , CaCO 0 3 ) on S3 buffering the local supersaturation. To do this, the initial speciation of the inlet solution was evaluated using the PHREEQC code and database, S3 and local equilibrium speciation was subsequently maintained throughout the simulated fluid cell. Due to the computational cost of this model, we were limited to a 2D plane parallel to the flow and perpendicular to the surface of the crystal. Initial conditions for the 2D model came from projections of slices of a 3D model that considered the transport of all species but without equilibrium reactions.
Good correspondence between the 2D and 3D models, especially close to the crystal surface, gave confidence in this approach. The effect of including equilibrium speciation within the 2D model had only a small (< 10%) effect on the final surface supersaturation, validating the limited species model adopted.

Mass Transport Model
Transport of component i was described by the stationary convection-diffusion equation where the flux J i is determined by the concentration c i and the diffusion coefficient D i , S2) and the solution velocity w is described by the stationary Navier-Stokes equation and the continuity equation where p is the pressure, ρ is the density, and µ = 8.925528 × 10 −4 Pa · s is the dynamic viscosity.
A no flux condition was applied to all boundaries except: 1. the inlet (the green boundary in Fig. S1a), where the inlet concentrations were fixed at the equilibrium values determined from the speciation model, and 2. the crystal surface (pink boundary in Fig. S1b), where a flux condition was instead applied (see the next section).
A no slip condition was applied to all boundaries, except: 1. the inlet, where the solution velocity corresponding to fully developed laminar flow at a given flow rate was specified, and 2. the outlet where the pressure was defined (p = 0).

Boundary Layer Thickness
The boundary layer thickness δ i associated with species i is defined such that the net flux of solute i from bulk (concentration c bulk,i ) to the surface (concentration c surf,i ) is D i (c bulk,i − c surf,i )/δ i . This net flux must balance with the crystal growth rate R/ω where R is the normal growth speed and ω is the molar volume of calcite, resulting in the equation where the average boundary layer δ = (δ Ca 2+ + δ CO 2− 3 )/2 is reported in the main text. To evaluate this parameter under realistic conditions, a flux boundary condition was applied to all faces of the crystal, with a normal flux proportional to k(S surf − 1) 2 for a value k chosen to approximately reproduce the normal growth speed reported experimentally. S1

S5
Note that this flux was evaluated locally, producing a non-uniform growth rate across the crystal surface-representative of the experimental crystals which featured an array of screw dislocations across the surfaces.
The boundary layer thickness in the finite element model exhibited a power law dependence δ = 418u −0.27 μm on the flow rate u (in mL/hr). This power law can reveal how the surface supersaturation in an AFM experiment will depend on the flow rate. Suppose the crystal surface structure has fully equilibrated to the solution, and that the normal growth rate of the surface takes the form R(S surf ) ≈ k(S surf − 1) 2 (where k = 1 nm/s for calcite, see Section 2.5). The experimental step velocities will be proportional to which follows from equation (S14), where ξ = kδγ ωDK 1/2 sp (S7) and γ is the activity coefficient.
In the AFM experiments, the step velocities were reported to plateau across a flow rate 30 u 40 mL/hr for nominal supersaturations as high as S bulk = 2. S4 Using the boundary layer thickness δ(u) obtained from the finite element model, it follows from equation (S6) that the step velocities will not vary by more than 3% across this flow range (even less at lower supersaturation), consistent with the change observed in the finite element model.
By contrast, AFM velocity measurements typically have errors of at least 10%, and suffer from stochastic history-dependence under mixed control if steady state is not fully realised, meaning that this flow rate dependence will be too weak to be resolved.

S6
Here we summarise our analysis of the AFM measurements of calcite reported in refs. S1,S5,S6

Critical lengths of nonequivalent step types
The geometry of calcite gives rise to two distinct step types, the acute and the obtuse step, which exhibit distinct critical lengths. S5 At first glance, it is surprising that they have different critical lengths. Each step type has the same pair of adjacent steps (acute on one side, obtuse on the other), and so the free energy cost φ in equation (1) should ostensibly be the same in both cases. The prevailing explanation for the different critical lengths, first presented in ref., S5 is that when a step segment emerges from a screw dislocation, a new corner has to be created. This corner may be curved (due to entropy), and the advancement of the step segment is limited by the cost of creating this curvature. The free energy cost of creating this curved corner will differ for the two step types, and so they will exhibit distinct critical lengths.
However, we believe this explanation is incorrect. A corner becomes curved because a curved corner is a lower free energy configuration than a sharp corner. Consequently, the advancement of a step segment will incur a greater free energy cost after the curved corner has fully emerged than beforehand. This latter cost will be responsible for the observed critical length, and it will be the same for both step types for the reason outlined above.
An alternative explanation for the different critical lengths comes from considering the different step type stabilities. In the absence of impurities, the velocity of one step type (of infinite length) will never equal zero under the exact same conditions as any other step type (of infinite length), as this would require a remarkable conspiracy between the elementary events associated with each step type. Consequently, as the supersaturation at a crystal/solution interface is decreased, one particular step type-the least stable-will be Growing Stationary Figure S2: (a) KMC model of a crystal with two nonequivalent step types (orange and blue). The only difference between the step types is the kink flux. (b) The dependence of (infinite) step velocity on S for parameters β = 2/7, γ = 1/2. (c) A simulation of an island confined to a channel of finite width (parameters β = 2/7, γ = 1/2, S = 1.2). The channel width is chosen to be the critical length of the orange step. The blue step nevertheless grows, demonstrating that the step types have distinct critical lengths.
first to halt. For calcite, the least stable step is the obtuse step. S4 We argue that the stability of a step type must partly determine the corresponding critical length.
To demonstrate this, we have performed a kinetic Monte Carlo simulation (details in Section 3.4) of a simple cubic crystal with two distinct step types, coloured orange and blue in Fig. S2(a). The elementary rate constants have been chosen so that each step type has the same step free energy, but the (infinite) step velocities have different dependencies on the supersaturation ( Fig. S2(b)). In particular, the orange step arrests at a higher supersaturation and is therefore less stable than the blue step. In our simulation, we confined an island S8 to a one-dimensional channel of fixed width. The width was chosen to be the critical length of the orange step. Over the course of the simulation, the orange step fluctuated about its starting position, but the blue step grew (Fig. S2(c)). It follows that the more stable (blue) step has a smaller critical length than the less stable (orange) step, and so the critical length of a step does depend on its stability.
The critical length given by equation (1) in the main text will correspond to the least stable (obtuse) step. To see this, consider a hypothetical crystal that is identical to calcite except the molecular kinetics at the acute step have been modified so that the acute step has the same stability as the obtuse step, but it is modified in a way that preserves its step free energy (this is possible since step free energy and step velocity are independent properties, see Section 3.4). Our claim then follows from three observations: 1. The obtuse step would have the same critical length in both calcite and the hypothetical crystal since the molecular kinetics for the obtuse step and the step free energy for the acute step are identical in both crystals.
2. The hypothetical crystal will have the same solubility as calcite since a crystal (calcite, in this case) ceases to grow once its least stable step (obtuse step) ceases to grow. Both crystals would therefore reach equilibrium under the same conditions.
3. In the hypothetical crystal, the two step types have the same stability (by definition), and so there is no doubt that equation (1) applies.
Note that point 2 in this argument is where the symmetry between the acute and the obtuse step gets broken.

Recovering L c Geometrically
Critical step lengths associated with a screw dislocation will be geometrically related to the terrace widths and step velocities of the accompanying step outcrop. This relationship will S9 depend on the crystal geometry.
In the case of calcite, and assuming a unit Burgers vector, the critical length of the obtuse step L c may be recovered geometrically from the step velocities (v ± ) and terrace widths (λ ± ) of the acute (-) and obtuse (+) steps. First, the average critical length L c of the acute and obtuse steps can be computed S5 via S9) measures the angle between adjacent spiral turns. S6 The obtuse critical length L c is then related to L c by the empirical relationship S5 The accuracy with which S bulk can be mapped to S surf using L c will depend on how accurately the step free energy φ is known. The step free energy of calcite has not been measured, but it can be estimated by either the surface free energy or the kink free energy. Measurements for the surface free energy include φ = 3k B T obtained from nucleation experiments, S7 and φ = 3.2k B T obtained using the vapour adsorption technique. S8 Estimates for the kink free energy can be derived from the observed kink density of calcite S9 using a model that relates the two; the answer depends on the model used, S10-S12 but falls between 2.5 φ/k B T 3. Bringing these various estimates together, we suggest that φ = (3 ± 0.5)k B T is a reasonable measure of the step free energy (in more familiar units, φ ≈ 120 ± 20 mJ/m 2 ). The shaded error bars in Fig. S3 show the degree of variation in our theoretical fits to the step velocities that arise from this uncertainty in step free energy.

S10
Significantly, the resulting error in step velocity (∼ 1 nm/s) is comparable in magnitude to the errors already present in the AFM velocity measurements.
Note that the data from which φ has been estimated are unlikely to suffer from the same mass transport distortions that we are ultimately addressing. In the nucleation experiments, S7 φ was calculated from nucleation rates which depend on the thermodynamics of critical clusters. However, critical-sized clusters are too small to be limited by diffusion.
This can be seen with a back-of-the-envelope calculation: In Section 2.4, we introduce a parameter α to characterise mass transport, and we estimate for calcite that α ∼ 0.1 for a boundary layer thickness δ ∼ 100 μm. It follows from equation (S14) that, for the same surface supersaturation, α will only deviate from 1 (which denotes surface control) when δ exceeds ∼ 1 μm. However, the boundary layer of a crystal in a stagnant solution equals the crystal radius, and micron-sized crystals are post-critical as they are observed to be stable in solution. The other experimental data were recorded at equilibrium.

Acute and Obtuse Step Velocities
The step velocities for the acute and obtuse steps, recalibrated against S surf , are both shown in Fig. S3. The solid line was fitted to the obtuse measurements via least squares regression subject to the constraint that the velocity must vanish at S surf = 1, yielding v + = 64.9(S surf − 1) nm/s (S11) while the dashed line was fitted to the acute measurements without constraint, v − = 16.5(S surf − 0.9) nm/s (S12) The lower stability of the acute step reported here is consistent with AFM observations. S4 Step velocity (nm/s) Obtuse v + Acute v Figure S3: The acute and obtuse step velocities of calcite, measured by AFM S1,S6 and recalibrated against S surf . Least squares regression lines are shown. The shaded error bars correspond to an uncertainty of ±0.5k B T in the step free energy φ. The acute step is evidently less stable than the obtuse step.

The Dependence of L c on S bulk
The dependence of the critical step length L c on the surface supersaturation S surf is straightforwardly described by equation (1) in the main text. However, the dependence of L c on S bulk under mixed surface/transport control will depend on how S surf and S bulk are related, which itself will depend on the prevailing distribution of steps across the surface. To explore this relationship, we consider two limiting cases.
First, suppose L c is measured against a series of supersaturations S bulk such that the surface is allowed to fully relax and achieve steady state with the solution for each individual measurement. Under this condition, the normal growth rate R(S surf ) will be the steady state function derived in Section 2.5. Furthermore, we approximate the Ca 2+ and CO 2− 3 ions as indistinguishable, with a shared diffusion coefficient D = (D Ca 2+ + D CO 2− 3 )/2 such that the boundary layer thickness (cf. equation (S5)) becomes

S12
where γ = 0.37 is the activity coefficient for the AFM experiments S1 computed from the Davies equation, and K sp is the solubility product of calcite.
If the boundary layer thickness δ is known, then S surf can be recovered for a given S bulk by numerically solving equation (S13). The critical length L c then follows from S surf via equation (1). The solid line in Fig. S4 shows the resulting steady state relationship between L c and S bulk for a boundary layer thickness δ = 108 μm, which comes from fitting to the three lowest-L c data points. We justify this fit on the basis that the L c measurements vary significantly across the S bulk range up until these three points where they appear to converge. This is significant since the surface should reconstruct most extensively for these three points since they correspond to the highest supersaturations. Moreover, the resulting functional form is consistent with our theoretical prediction for steady state, and the optimised δ value is in reasonable agreement with the δ ≈ 160 μm obtained from our finite element model.
Note that all points that deviate from this line have failed to achieve steady state and are therefore sensitive to their surface histories. Longer equilibration times would have moved each point nearer to the steady state curve.
In the second case, suppose L c is measured across a range of S bulk , but over a sufficiently short time-frame that the distribution of steps across the surface remains approximately unchanged across all measurements. The significance of not allowing the surface to reconstruct is that it breaks the feedback between step density and surface supersaturation, yielding a characteristic relationship between L c and S bulk . To analyse this case, it helps to introduce the parameter which follows from equation (S13). This parameter characterises mass transport, with α = 1 denoting surface control and a smaller α denoting greater transport control. For a static surface structure, the normal growth rate is linearly proportional to the relative supersaturation, R(S surf ) ∼ (S surf − 1), and so α will be constant with a value determined by the history-dependent surface structure (cf. equation (S14)).
Utilising α, the critical length can be shown to simplify to a linear function of (ln S bulk ) −1 at low supersaturation, where equation (S17) employs a Taylor expansion, and equation (S18) follows from the S14 binomial approximation. Note that equation (S18) is exact when α = 1, as well as in the limit S bulk → 1+ irrespective of α. The largest error occurs when α 1 and S bulk is large.
The largest supersaturation sampled in the AFM measurements was S bulk = 2 for which the approximation (equation (S18)) underestimates the critical length by 6% at worst.
The critical length has a larger gradient (α −1 × larger) when plotted against (ln S bulk ) −1 (equation (S18)) than when plotted against (ln S surf ) −1 (equation (S15)). Moreover, when L c is plotted against (ln S bulk ) −1 it is offset from the origin: L c is projected to vanish at ln(S bulk ) −1 = 1 2 (1 − α) (≈ 1 2 since α 1 under mixed control). Optimising α = 0.0817 to fit equation (S18) to the measurements from ref. S5 reproduces both the gradient and the offset (compare the dashed line to the open circles in Fig. S4), validating this analysis.

Steady State Normal Growth Rate
The normal growth rate of calcite R(S surf ) under steady state conditions (where the surface structure has fully equilibrated to the prevailing solution) can be computed as follows.
2. Compute the obtuse critical length L c from equation (1).
3. Compute the average critical length L c from equation (S10).
5. The normal growth rate is then R(S surf ) = av ± /λ ± The resulting function is plotted in Fig. 2c of the main text (solid line). Note that the growth rate that results from this procedure can be approximated by R(S surf ) ≈ k(S surf − 1) 2 , k = 1 nm/s (S19) S15 where one (S surf − 1) term comes from the step velocity, and a second comes from the step density via the approximation ln S surf ≈ S surf − 1.

Scale of Solution Inhomogeneities
Under mixed surface/transport control, a crystal surface with a spatially non-uniform growth rate will induce a spatially non-uniform supersaturation across the surface of the crystal.
Acting against these gradients, however, is the lateral diffusion of solutes. The magnitude of solution inhomogeneity across a crystal surface will therefore depend on the balance between the rate of crystal growth and the rate of solute diffusion.
As a crude estimate of this balance, the time-scale over which the steps advance by a monomolecular row (of length a) is τ = a/v, and in that window of time, solutes will diffuse a distance l ∼ √ 4Dτ laterally across the surface. Using a = 0.32 nm, v ∼ 10 nm/s, and D ∼ 10 −9 m 2 /s reveals l ∼ 10 µm for calcite. In other words, the solution locally mixes over a scale of ∼ 10 µm.
In support of this number, Bracco et al. S13 reported that calcite step velocity measurements required 20 minutes to equilibrate, during which the steps travelled at an average speed of approximately 10 nm/s. This equates to ∼ 10 µm of surface reconstruction before the surface (in the vicinity of a screw dislocation) lost its sense of history.

Kinetic Monte Carlo Simulation
To study the dependence of step velocity on step length, we developed a simple Kossel model with the same step free energy as calcite. The model was integrated through time using KMC simulation.

S16
S surf 1 Figure S5: Reaction rates in the Kossel model of calcite. The attachment (green) and detachment (red) rates depend on the number of bonds (shaded squares). Reactions for zero-and four-bonded sites are described in the text.

Kossel Model of Calcite
In our model for a propagating crystal step, the crystal is represented by a two-dimensional grid. Each site is either occupied by an ion or it is not, and has between 0 and 4 (inclusive) neighbouring bonds.
Each site will transition between occupied and unoccupied states at a rate determined by the number of bonds it has, as summarised in Fig. S5. In general, each individual reaction has been constructed so that (i) the creation of a unit of surface induces a free energy cost φ = 3k B T , and (ii) solute attachment reduces the free energy by ∆µ = k B T ln S surf .
Enforcing the thermodynamic definitions of S surf and φ at the level of each individual reaction guarantees that they also hold at the level of the entire step.
As an approximation, we omit this thermodynamic equilibrium for sites with zero and four bonds. In our model, solutes dissolve from zero-neighbour sites at a rate 1 and attach to four-neighbour sites at a rate S. In both cases, the reverse reaction should occur at a rate e −4φ/k B T to maintain the desired equilibrium. However, since these reverse reactions are rare and occur away from the step edge, we exclude them with the approximation e −4φ/k B T ≈ 0.

Measuring Step Velocities
The crystal configuration is represented on a two-dimensional boolean-valued grid of dimension w × h. The height is h = 2 12 in all cases, and is terminated at both ends by reflective S17 boundary conditions. When measuring a step segment of length L, the grid width is set to w = L/a with aperiodic boundaries. When estimating the velocity of an (infinite) step, the grid width is set to w = 2 10 with periodic boundaries.
The crystal initially occupies the area [0, w) × [0, h/2) so that the reflective boundaries ensure there is a single growth front. The configuration is integrated over time using KMC. S14 The advancing crystal is periodically displaced back to the centre of the grid (along y) to prevent it from exhausting the available grid space.
The system is equilibrated for 10 9 integration steps, followed by a production period of 10 9 steps over which the velocity is measured. This routine is repeated 100 times, yielding an average velocity with an aggregate simulation time of 10 11 integration steps, plus a measure of error. The largest standard error associated with any normalised velocity measurement was ±0.18%. The resulting KMC measurements (plotted in Fig. 3 of the main text) are normalised with respective to the theoretical critical length L c (equation (1) in the main text). Note that the step velocity vanishes at the theoretical critical length, as expected.

Replication of AFM
Teng et al. S5 measured the dependence of step velocity on step length, v(L), across a nominal supersaturation range 1.19 ≤ S bulk ≤ 1.51 using AFM. The resulting velocity profile exhibited no discernible dependence on supersaturation across this range, and is shown in Fig. 3 of the main text.
We computed v(L) using our KMC model for two cases: 1. surface control, where S bulk = S surf = 1.51, and 2. mixed surface/transport control, where S bulk = 1.51 but S surf = 1.043 in accordance with the experimental critical length L c ≈ 45.9 nm that corresponded to this nominal supersaturation (cf. equation (1)). S18

Model of nonequivalent step types
In Section 2.1, a KMC model of nonequivalent step types was briefly outlined. We elaborate on it here.
The model is a modified version of that previously described (Section 3.1, Fig. S5).
Instead of all four step orientations of an island being identical, a pair of steps (blue in Fig. S2) are given distinct elementary rate constants. For both step types (blue and orange), the creation and annihilation of kinks occurs at identical rates, and so they share a step free energy. However, while solutes attach to (detach from) orange kinks at a rate S (1), they attach to (detach from) blue kinks at a rate βS (βγ). Consequently, the orange step velocity will scale as ∼ (S − 1) while the blue step velocity will approximately scale as ∼ β(S − γ). This allows both the magnitude of the step velocity and the stability of the step to be adjusted. For the simulation in Fig. S2(c), we chose β = 2/7, γ = 1/2, S = 1.2 so that the two step types have the same (infinite) step velocity, demonstrating that their different velocities within the channel result from their distinct length-dependencies, i.e. their distinct critical lengths.